#Authors: Natili, Pellegata, Visconti
#Date: 06 November 2025
#Data: Eco-Social UNIMI (Ferrera et al. 2022)
#Paper Title: What Kind of Energy Transition? Public Opinion Trade-Offs among Economic
#       Growth, Ecological Sustainability and Equity
#Journal: European Journal of Political research
#Object code: Appendix figures

#------------------------------------------------------------------------------#
#LOAD/INSTALL REQUIRED PACKAGES
#------------------------------------------------------------------------------#

#Load the relevant packages or install them if not available
if (!require("pacman")) install.packages("pacman")
pacman::p_load(cregg, tidyverse, ggplot2, patchwork)

#If cregg does not install please use the the following code:
#remotes::install_github("leeper/cregg")

#------------------------------------------------------------------------------#
#SPECIFY PATH TO WORKING DIRECTORY WITH DATASET
#------------------------------------------------------------------------------#
#Set working directory Windows
#setwd("C:\Users\fvisconti\Documents\Replication")

#Set working directory Mac OS
setwd("/Users/fv/Library/CloudStorage/Dropbox/ActiveWorks/myREScUE/Papers/Eco-social/eco-social-shared/Paper/EJPR-submission/R2/Replication")
#------------------------------------------------------------------------------#

#Clear all objects in environment
rm(list = ls())

#------------------------------------------------------------------------------#
#LOAD DATASET
#------------------------------------------------------------------------------#
#Load data object in long format ready for analysis
load(file = "dataset.RData")

#------------------------------------------------------------------------------#
#CREATE APPENDIX FIGURES FOLDER
#------------------------------------------------------------------------------#

dir.create("appendix_figures", showWarnings = FALSE)


#------------------------------------------------------------------------------#
##### FIGURE A1 - COUNTRY RESULTS#####
#------------------------------------------------------------------------------#

cj_anova(subset(long, !is.na(long$qcountry)), chosen_rec ~ Energy*Compensation + Financing, id = ~caseid, by = ~qcountry )
cj_anova(subset(long, !is.na(long$qcountry)), chosen_rec ~ Energy, id = ~caseid, by = ~qcountry )
cj_anova(subset(long, !is.na(long$qcountry)), chosen_rec ~ Compensation, id = ~caseid, by = ~qcountry )
cj_anova(subset(long, !is.na(long$qcountry)), chosen_rec ~ Financing, id = ~caseid, by = ~qcountry )

#Compute marginal means
cnt1 <- cj(
  subset(long, !is.na(long$qcountry)),
  chosen_rec ~ Energy*Compensation + Financing, id = ~caseid, by = ~ qcountry,
  weights = ~ weight, estimate = "mm"
)

#Define plot text style
myFaces <- c(rep('plain', 4), "bold",
             rep('plain', 4), "bold",
             rep('plain', 3), "bold")

#Define plot
cnt1plot <- plot(cnt1, group = "qcountry",
                 vline = 0.5, 
                 header_fmt = "%s",
                 size = 3) +
  aes(shape = feature, group = "feature") + # map group to `shape` aesthetic
  scale_shape_manual(values=c(15, 16, 17)) +
  labs(title = "Support for policy package by country") +
  guides(colour=guide_legend(title="Country")) +
  ggplot2::theme(plot.title = element_text(size=20),
                 plot.subtitle = element_text(size=18),
                 axis.text=element_text(size=14),
                 axis.title=element_text(size=16),
                 axis.text.y = element_text(face=myFaces),
                 text = element_text(size = 18), 
                 legend.position = "none") +
  facet_wrap(~BY, ncol = 4L) +
  ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower",
                                              xmax = "upper"),
                          linewidth = 0.75, height = 0, na.rm = TRUE)

cnt1plot

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A1.png',
       plot = cnt1plot, width = 20, height = 12)

#------------------------------------------------------------------------------#
##### FIGURE A2 - SUPPORT FOR VARIATIONS OF JUST TRANSITION PACKAGES #####
#------------------------------------------------------------------------------#

#Create interaction between energy and compensations variables
long$Interaction1 <- interaction(long$Energy, long$Compensation)

#Compute marginal means for all combinations of packages
interfit3 <- cj(long, chosen_rec ~ Interaction1 ,
                id = ~ caseid,
                weights = ~ weight,
                estimate = "mm",
                h0 = 0.5,
                by = ~Financing)

#JUST TRANSITION VARIATIONS
jtdat <- interfit3 %>%
  filter(
    (BY == "Wealth tax" & level == "Invest in renewables.Support affected workers") | 
      (BY == "Carbon tax" & level == "Invest in renewables.Support affected workers") | 
      (BY == "Wealth tax" & level == "Invest in renewables.Cash to low-income HHs")  |
      (BY == "Carbon tax" & level == "Invest in renewables.Cash to low-income HHs") 
  )

#Assign name to policy packages
jtdat$name <- NA

jtdat$name[jtdat$BY == "Wealth tax" & 
             jtdat$level == "Invest in renewables.Support affected workers"] <- 
  "JUST TRANSITION\nInvest in renewables\nSupport workers\nWealth tax"

jtdat$name[jtdat$BY == "Carbon tax" & 
             jtdat$level == "Invest in renewables.Support affected workers"] <-
  "JUST TRANSITION\nInvest in renewables\nSupport workers\nCarbon tax"

jtdat$name[jtdat$BY == "Wealth tax" & 
             jtdat$level == "Invest in renewables.Cash to low-income HHs"] <- 
  "JUST TRANSITION\nInvest in renewables\nCash to low-income HHs\nWealth tax"

jtdat$name[jtdat$BY == "Carbon tax" & 
             jtdat$level == "Invest in renewables.Cash to low-income HHs"] <-
  "JUST TRANSITION\nInvest in renewables\nCash to low-income HHs\nCarbon tax"

#Transform "name" to factor variable
jtdat$name <- factor(jtdat$name,
                     levels = c("JUST TRANSITION\nInvest in renewables\nSupport workers\nWealth tax",
                                "JUST TRANSITION\nInvest in renewables\nSupport workers\nCarbon tax",
                                "JUST TRANSITION\nInvest in renewables\nCash to low-income HHs\nWealth tax",
                                "JUST TRANSITION\nInvest in renewables\nCash to low-income HHs\nCarbon tax"))

table(jtdat$name)

#Reorder levels based on values
jtdat <- jtdat %>%
  mutate(name = reorder(name, -estimate)) # Reorder name based on estimate

#Define the plot in Black and White
jtplot <- ggplot(jtdat, aes(x = name, y = estimate)) +
  geom_bar(stat = "identity", position = position_dodge(), colour = "black", fill = NA) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, position = position_dodge(0.9)) +
  labs(x = "Policy paradigm", y = "Marginal mean estimate") +
  theme_classic() +
  theme(legend.position = "none",
        text = element_text(size = 16), 
        axis.text = element_text(size = 12)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black")

jtplot

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A2.png',
       plot = jtplot, width = 10, height = 6)

#------------------------------------------------------------------------------#
##### FIGURE A3 - SUPPORT FOR VARIATIONS OF POST-GROWTH PACKAGES #####
#------------------------------------------------------------------------------#

#POST GROWTH VARIATIONS
pgdat <- interfit3 %>%
  filter(
    (BY == "Wealth tax" & level == "Downscale production.Cash to all HHs") | 
      (BY == "Increasing public debt" & level == "Downscale production.Cash to all HHs") |
      (BY == "Carbon tax" & level == "Downscale production.Cash to all HHs")
  )

#Assign name to policy packages
pgdat$name <- NA

pgdat$name[pgdat$BY == "Wealth tax" & 
             pgdat$level == "Downscale production.Cash to all HHs"] <- "POST GROWTH\nDownscale production\nCash all HHs\nWealth tax"

pgdat$name[pgdat$BY == "Increasing public debt" & 
             pgdat$level == "Downscale production.Cash to all HHs"] <- "POST GROWTH\nDownscale production\nCash all HHs\nPublic debt"

pgdat$name[pgdat$BY == "Carbon tax" & 
             pgdat$level == "Downscale production.Cash to all HHs"] <- "POST GROWTH\nDownscale production\nCash all HHs\nCarbon tax"

#Transform "name" to factor variable
pgdat$name <- factor(pgdat$name,
                     levels = c("POST GROWTH\nDownscale production\nCash all HHs\nWealth tax",
                                "POST GROWTH\nDownscale production\nCash all HHs\nPublic debt",
                                "POST GROWTH\nDownscale production\nCash all HHs\nCarbon tax"))

table(pgdat$name)

#Reorder levels based on values
pgdat <- pgdat %>%
  mutate(name = reorder(name, -estimate)) # Reorder name based on estimate

#define the plot in Black and White
pgplot <- ggplot(pgdat, aes(x = name, y = estimate)) +
  geom_bar(stat = "identity", position = position_dodge(), colour = "black", fill = NA) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, position = position_dodge(0.9)) +
  labs(x = "Policy paradigm", y = "Marginal mean estimate") +
  theme_classic() +
  theme(legend.position = "none",
        text = element_text(size = 16), 
        axis.text = element_text(size = 12)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black")

pgplot

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A3.png',
       plot = pgplot, width = 8, height = 8)

#------------------------------------------------------------------------------#
##### FIGURE A4 - SUPPORT FOR POLICY PACKAGES BY OCCUPATION #####
#------------------------------------------------------------------------------#

# Identify all unique values of the occ variable
unique_occ <- unique(long$occ2)

# Initialize an empty data frame to store the combined dataset
combined_data <- data.frame()

# Loop over each unique value to create and store plots
for(occ2 in unique_occ) {
  
  # Subset the data for the current country
  subset_long <- long[long$occ2 == occ2,]
  
  # Perform your analysis here (cj, filter, name assignment, etc.) similar to the provided code
  # Adapted for subset_long instead of longFra
  interfit4 <- cj(subset_long, chosen_rec ~ Interaction1 ,
                  id = ~ caseid,
                  weights = ~ weight,
                  estimate = "mm",
                  h0 = 0.5,
                  by = ~Financing)
  
  
  #KEEP POLICY PARADIGMS SPECIFIED IN PAPER
  packagedat4 <- interfit4 %>%
    filter(
      (BY == "Cuts in social spending" & level == "Rely on fossil fuels.Support affected companies") | #GROWTH FIRST
        (BY == "Cuts in social spending" & level == "Invest in renewables.Support affected companies") | #GREEN GROWTH
        (BY == "Wealth tax" & (level == "Invest in renewables.Support affected workers")) | #JUST TRANSITION
        (BY == "Wealth tax" & level == "Downscale production.Cash to all HHs") #POST GROWTH
    )
  
  #Assign name to policy packages
  packagedat4$name <- NA
  
  
  packagedat4$name[packagedat4$BY == "Cuts in social spending" & 
                     packagedat4$level == "Rely on fossil fuels.Support affected companies"] <- "GROWTH FIRST\nFossil fuels\nSupport companies\nSocial cuts"
  
  
  packagedat4$name[packagedat4$BY == "Cuts in social spending" & 
                     packagedat4$level == "Invest in renewables.Support affected companies"] <- "GREEN GROWTH\nInvest in renewables\nSupport companies\nSocial cuts"
  
  
  packagedat4$name[packagedat4$BY == "Wealth tax" & 
                     packagedat4$level == "Invest in renewables.Support affected workers"] <- "JUST TRANSITION\nInvest in renewables\nSupport workers\nWealth tax"
  
  packagedat4$name[packagedat4$BY == "Wealth tax" & 
                     packagedat4$level == "Downscale production.Cash to all HHs"] <- "POST GROWTH\nDownscale production\nCash all HHs\nWealth tax"
  
  
  #Transform "name" to factor variable
  packagedat4$name <- factor(packagedat4$name,
                             levels = c("GROWTH FIRST\nFossil fuels\nSupport companies\nSocial cuts",
                                        "GREEN GROWTH\nInvest in renewables\nSupport companies\nSocial cuts",
                                        "JUST TRANSITION\nInvest in renewables\nSupport workers\nWealth tax",
                                        "POST GROWTH\nDownscale production\nCash all HHs\nWealth tax"))
  
  table(packagedat4$name)
  
  packagedat4$occ2 <- occ2  # Add occupation information
  
  combined_data <- rbind(combined_data, packagedat4)
  
}


#Exclude Not in work/Others and Homemaker
combined_data <- combined_data %>%
  filter(!occ2 %in% c("Homemaker", "Not in work/Others", "Student"))

#Define the plot
occplot <- ggplot(combined_data, aes(x = name, y = estimate)) +
  geom_bar(stat = "identity", position = position_dodge(), colour = "black", fill = NA) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, position = position_dodge(0.9)) +
  labs(x = "Policy paradigm", y = "Marginal mean estimate") +
  facet_wrap(.~occ2) +
  theme_classic() +
  theme(legend.position = "none",
        text = element_text(size = 12), 
        axis.text = element_text(size = 8),
        #axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1)
  ) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black")

occplot

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A4.png',
       plot = occplot, width = 14, height = 10)

#------------------------------------------------------------------------------#
##### Diagnostics & Robustness checks#####
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
#FIGURE A5 - Frequencies of conjoint features
#------------------------------------------------------------------------------#

p <- plot(cj_freqs(long, chosen_rec ~ Energy*Compensation + Financing , id = ~caseid))
p

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A5.png',
       plot = p, width = 6, height = 8)

#------------------------------------------------------------------------------#
#FIGURE A6 - #Differences left-right (Profile A Vs B)
#------------------------------------------------------------------------------#

#Define interaction
f1 <- chosen_rec ~ Energy*Compensation + Financing

long$profile_fac <- factor(long$profile)

#Define the plot
p <- plot(cj(long, f1, id = ~caseid, by = ~profile_fac, estimate = "mm"),
          group = "profile_fac", vline = 0.5,
          weights = ~ weight, size = 3) +
  theme(plot.title = element_text(size=20),
        plot.subtitle = element_text(size=16),
        legend.position = "bottom",
        legend.title = element_blank(), 
        legend.text=element_text(size=16),
        axis.text=element_text(size=16),
        axis.title=element_text(size=16)) +
  labs(title = "Conjoint experiment: effects")
p

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A6.png',
       plot = p, width = 10, height = 10)

#------------------------------------------------------------------------------#
#FIGURE A7 - #Differences Task 1 and Task 2
#------------------------------------------------------------------------------#

#Define the factor for tasks for computing marginal means
long$task <- factor(long$task,
                    levels = c(1, 2),
                    labels = c("Task 1", "Task2"))

#Compute marginal means
fittask <- cj(long, f1, id = ~ caseid, estimate = "mm", by = ~task,
              weights = ~weight)

#Define the plot
p <- plot(fittask, id = ~caseid, group = "task", vline = 0.5, size = 3) +
  theme(plot.title = element_text(size=20),
        plot.subtitle = element_text(size=16),
        legend.position = "bottom",
        legend.title = element_blank(), 
        legend.text=element_text(size=16),
        axis.text=element_text(size=16),
        axis.title=element_text(size=16)) +
  labs(title = "Conjoint experiment: effects")

p

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A7.png',
       plot = p, width = 10, height = 10)


#------------------------------------------------------------------------------#
#FIGURE A8 - #Marginal means computed on rating variables
#------------------------------------------------------------------------------#

#Create rating 2 variable without missing values
long$rating2 <- long$rating

#Drop NAs from rating
table(long$rating2)

long$rating2[long$rating2 == 977] <- NA

#Rescale rating 2 0-1
long$rating2 <- scales::rescale(long$rating2)

summary(long$rating2)

#Compute the marginal means
fitrating <- cj(long, rating2 ~ Energy*Compensation + Financing, id = ~ caseid,
                weights = ~ weight, estimate = "mm")

#Define the plot
p <- plot(fitrating, vline = .5) +
  theme(plot.title = element_text(size=20),
        plot.subtitle = element_text(size=16),
        legend.position = "bottom",
        legend.title = element_blank(), 
        legend.text=element_text(size=16),
        axis.text=element_text(size=16),
        axis.title=element_text(size=16)) +
  geom_point(size=4) +
  labs(title = "Conjoint experiment: effects"
  )

p

#Save the plot
ggsave(filename = 'appendix_figures/Figure_A8.png',
       plot = p, width = 10, height = 10)

#------------------------------------------------------------------------------#